Navier-Stokes hydrodynamics of thermal collapse in a freely cooling granular gas 
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We employ Navier-Stokes granular hydrodynamics to investigate the long-time behavior of clus- 
tering instability in a freely cooling dilute granular gas in two dimensions. We find that, in circular 
containers, the homogeneous cooling state (HCS) of the gas loses its stability via a sub-critical pitch- 
fork bifurcation. There are no time-independent solutions for the gas density in the supercritical 
region, and we present analytical and numerical evidence that the gas develops thermal collapse 
unarrested by heat diffusion. To get more insight, we switch to a simpler geometry of a narrow- 
sector-shaped container. Here the HCS loses its stability via a transcritical bifurcation. For some 
initial conditions a time-independent inhomogeneous density profile sets in, qualitatively similar to 
that previously found in a narrow-channel geometry. For other initial conditions, however, the dilute 
gas develops thermal collapse unarrested by heat diffusion. We determine the dynamic scalings of 
the flow close to collapse analytically and verify them in hydrodynamic simulations. The results 
of this work imply that, in dimension higher than one, Navier-Stokes hydrodynamics of a dilute 
granular gas is prone to finite-time density blowups. This provides a natural explanation to the 
formation of densely packed clusters of particles in a variety of initially dilute granular flows. 

PACS numbers: 45.70.Qj, 47.20.Ky 



I. INTRODUCTION 

Finite-time singularities are ubiquitous in hydrody- 
namic flows, and whenever they appear they attract 
much interest [1-3] . Wave breaking in ideal gas flows [4] , 
pinching of droplets [5] and creation of black holes [6] are 
examples of finite-time singularities which dramatically 
change the dynamics in the system when they appear. 
Singularities often signal breakdown of theory. In this 
case cither the theory is regularized by introducing ad- 
ditional macroscopic mechanisms, or the continuum de- 
scription is abandoned in favor of a discrete microscopic 
theory. A long-standing problem of hydrodynamics is 
whether a finite-time vorticity singularity develops in an 
initially smooth three-dimensional incompressible flow. 
Despite considerable effort, both the inviscid Euler case 
[7], and the viscous, Navier-Stokes case [8] remain unre- 
solved. 

The present paper deals with a somewhat simpler prob- 
lem of finite-time density singularities in dilute granular 
gases. In contrast to ordinary molecular gases, granu- 
lar gases "cool" spontaneously because of the inelastic 
collisions between the grains. Strikingly, even a slight in- 
elasticity of the collisions may cause the granular gas to 
form dense clusters of particles, even when starting from 
a dilute and macroscopically homogeneous initial state. 
This fascinating example of structure formation has been 
investigated by means of molecular dynamics simulations 
[9-24] and hydrodynamic simulations [20, 21, 23-30]. In 
addition, as many as five different hydrodynamical sce- 
narios have been suggested for the clustering instability; 
see Ref. [24] for a concise critical review. The cluster- 
ing instability in granular gases has an interesting analog 
in the form of condensation instability in optically thin 
gases and plasmas that cool by their own radiation [31]. 

The Navier-Stokes granular hydrodynamics is a nat- 



ural starting point for a physicist who wants to under- 
stand and predict collective physical phenomena in gran- 
ular gases [32, 33]. Granular hydrodynamics provides 
powerful insights into the physics of granular flow, not 
in a small part because its predictions are occasionally 
successful beyond its (quite restrictive) formal limits of 
validity [33]. In the present work we will employ the 
Navier-Stokes granular hydrodynamics in an attempt to 
better understand the physics of clustering, including its 
ultimate stage of thermal collapse, in a freely-cooling di- 
lute granular gas in two dimensions (2D). We will as- 
sume throughout the paper that the particle collisions 
are nearly elastic, the local gas density (that we denote 
by p) is much smaller than the close-packing density, and 
the Knudsen number of the granular gas is very small: 
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<Cl, pa d <l, and l free /L<^l. (1) 



Here < a < 1 is the coefficient of normal restitution 
of binary collisions, a is the particle diameter, d > 1 is 
the dimension of space, I free is the mean free path of 
the gas, and L is the characteristic hydrodynamic length 
scale. Under these assumptions (the second and third 
of them need to be verified a posteriori, after the hy- 
drodynamic problem in question is solved) the Navier- 
Stokes hydrodynamics provides a quantitatively accurate 
leading-order theory [32, 33]. 

In a sufficiently small container, a freely cooling gran- 
ular gas becomes homogeneous because of heat diffusion. 
From then on, the gas temperature decays according to 
Haff's law T ~ (to + t)~ 2 , where t is time, and = const 
[34]. This Homogeneous Cooling State (HCS) is linearly 
unstable, however, when the container size exceeds a crit- 
ical value [10, 13, 14, 35]. The instability manifests itself 
by the appearance of density inhomogeneities (the clus- 
tering instability) and vortical motions (the shear insta- 
bility) , each of them with its own critical system size. 
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Preceding our work is a series of papers [23, 24, 28- 
30, 36] which dealt with nonlinear clustering dynamics of 
a freely-cooling granular gas in a narrow channel. Fouxon 
et al. [29, 30] investigated the nonlinear clustering an- 
alytically and numerically, using ideal (Euler) granular 
hydrodynamics for a dilute gas which neglects heat dif- 
fusion and viscosity. They found a family of exact an- 
alytic solutions for the clustering flow. These solutions 
describe the evolution of an initially smooth flow toward 
thermal collapse: a finite-time density blowup. Close to 
the blowup time t r , the maximum gas density exhibits 
a power law behavior ~ (t* — t)~ 2 . The minimum tem- 
perature of the gas vanishes as (t* — t) 2 . The velocity 
gradient blows up as <~ (i* — whereas the veloc- 

ity itself remains continuous and develops a cusp, rather 
than a shock discontinuity, at the singularity. Molecular 
dynamics simulations of a freely cooling granular gas in 
a channel geometry [24] showed an excellent agreement 
with one of the exact solutions until sufficiently close to 
the attempted collapse, where heat diffusion becomes im- 
portant. At finite gravity, a one-dimensional (ID) cooling 
flow also develops thermal collapse, and exhibits the same 
asymptotic scaling behavior near the singularity, even in 
the framework of non-ideal, Navier-Stokes granular hy- 
drodynamics [23]. 

Despite the blowup of the density and the velocity gra- 
dient, the gas pressure behaves regularly in all exam- 
ples of thermal collapse considered in Refs. [23, 29, 30]. 
At zero gravity the collapse occurs at almost constant 
pressure [29, 30], whereas at finite gravity the pressure 
is approximately hydrostatic [23]. The authors of Refs. 
[23, 36] took advantage of these salient features to greatly 
simplify the full ID hydrodynamic problem, which takes 
a full account of heat diffusion and viscosity. It turns out 
that, at zero gravity, the heat diffusion arrests the ther- 
mal collapse in ID [36]. Instead of blowing up, the gas 
density approaches a time-independent Inhomogeneous 
Cooling Sate in which the heat diffusion balances the 
clustering instability brought about by the collisional en- 
ergy loss. 

These ID hydrodynamic results, however, are strik- 
ingly different from those observed in molecular dynam- 
ics simulations of a freely cooling granular gas in 2D 
[10, 11, 14, 16], where the gas density continues to grow 
until close-packed clusters of of particles are formed. To 
address this issue, the nonlinear hydrodynamic theory of 
clustering must be extended to higher dimensions. The 
first step in this direction has been recently taken by 
Fouxon [37]. Assuming a zero pressure gradient and us- 
ing Lagrangian coordinates he showed that, in the frame- 
work of ideal hydrodynamics, thermal collapse persists in 
any dimension. In Lagrangian coordinates the flow sin- 
gularity turns out to be identical to the one found in ID 
[36]. Fouxon [37] also briefly considered the non-ideal, 
Navier-Stokes flow. Here he found inhomogeneous so- 
lutions with time-independent density profiles, but left 
open the question of their stability. 

In this work we apply Navier-Stokes hydrodynamics to 



a freely-cooling granular gas in 2D, in a circular geome- 
try. The final outcome of our work is quite different both 
from the previous ID results [36], and from the results of 
Ref. [37]. The hydrodynamic problem is formulated in 
section II. In Section III we perform a marginal stability 
analysis of the HCS in a circular container and find the 
critical radii of the container for the clustering modes and 
shear modes. One of the surprising results, presented in 
section IV, shows that the clustering instability appears 
in this geometry via a sub-critical bifurcation. This sug- 
gests that the dilute flow can develop thermal collapse, 
and we support this conjecture by hydrodynamic simula- 
tion (that is, by solving the hydrodynamic equations nu- 
merically) . To get a more analytical insight we switch, in 
section V, to a narrow-sector-shaped container, see Fig. 
1. Here any flow structure in the azimuthal direction is 
suppressed by heat diffusion. As a result, the sector ge- 
ometry renders the simplicity of an effectively ID (purely 
radial) time-dependent flow, while retaining some critical 
2D features. We find that, in the sector geometry, the 
bifurcation of the clustering instability is transcritical. 
Here stable time-independent azimuthally-symmetric in- 
homogeneous density states exist with a density peak at 
the circumference of the sector. For these states the den- 
sity growth is exactly balanced by heat diffusion. But on 
the same bifurcation diagram there also regions which 
do not correspond to any time-independent state. Our 
hydrodynamic simulations in these regions give a strong 
evidence for thermal collapse unarrested by heat diffu- 
sion, this time at the vertex of the sector. We investigate 
thermal collapse analytically in section VI. We determine 
scaling laws, in time and in space, for an ideal collapse 
which we assume to develop on the background of an al- 
most constant pressure. Using the scaling laws we find 
that the non-ideal transport mechanisms - heat diffusion 
and viscosity - are in general unable to arrest the collapse 
in any dimension higher than one. 




FIG. 1: Two-dimensional container geometries considered in 
this work. 



II. HYDRODYNAMIC FORMULATION 

A. Governing equations and boundary conditions 

The granular hydrodynamic equations express local 
conservation of mass and momentum, and local energy 
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balance, in the granular gas. Assuming the strong in- 
equalities (1), one can use the following set of hydrody- 
namic equations [32, 33, 38]. The continuity equation, 



dp 
dt 



+ V • Gov) = 0, 



(2) 



and the momentum equation, 
dv 



dt 



+ (v • V)v 



Vp + ^ V-(VTS). (3) 



are the usual Navier- Stokes equations for a compressible 
flow. Here p(r,t) is the gas density, v(r, t) is the gas 
velocity, and ^oVTS is the viscosity tensor. In Cartesian 
coordinates 



Sfei = d k Vi + diV k - 5 ik V ■ v, 
see Re f. [4]. The energy equation reads: 



(4) 



dT 



— + (v V)T = -( 7 - 1)TV • v + -jV . (VTVT) 



Apl 



T 3/2 + M"f - ]WT ^ 2 



2p 



(5) 



where 7 is the adiabatic index of the gas (7 = 2 and 
5/3 for d — 2 and d = 3, respectively), and the ideal- 
gas equation of state p = pT is used. Equations (2)- (5) 
differ from the hydrodynamic equations for a dilute gas 
of elastically colliding spheres only by the presence of the 
inelastic energy loss term — ApT 3 / 2 which is proportional 
to the average energy loss per collision, ~ (1 — cn 2 )T, and 
to the collision rate, <~ pT 1 / 2 . In the viscous heating 
term of Eq. (5), X 2 is a scalar. In Cartesian coordinates 
it is the sum of the squared components (Eit) 2 . 

The transport coefficients in Eqs. (3) and (5) arc 
known from kinetic theory: A = 27r^ d_1 ^ 2 (l — 
a 2 )<T d - 1 /[dT(d/2)} (see e.g. [38]), where T(. . . ) is the 
gamma-function, and d > 2 is the dimension of space, 
so that d = 2 corresponds to disks, and d — 3 to hard 
spheres. Furthermore, v a = (2<Ty/n)^ 1 and n = 4^ in 
2D, and v = 5(3a 2 y/F)- 1 and k = 15u /8 in 3D [32]. 
Throughout this paper we assume that the particles have 
a unit mass. 

It is often convenient to rewrite the energy equation as 
an evolution equation for the pressure: 



dp 

+ (v • V)p = -7pV • v + k V 




A^/y/ 2 



p \p 

2 



(6) 



In a circular geometry that we will be dealing with 
throughout this paper, the natural coordinates are the 
polar coordinates r, 9. The tensor S is given in these 
coordinates by 



dv r 
dr 



v r 
r 



1 dve 
r d9 



dve 
dr 



1 dv r 
7~dT 



•(7) 



As the viscosity tensor is symmetric and traceless [4], 
the other components can be readily found: Yigg = 

The hydrodynamic equations need to be supplemented 
by boundary conditions. Let R be the radius of the con- 
tainer. We choose the boundary conditions so as not 
to introduce any sources or sinks of mass, momentum 
and energy, and also not to break azimuthal symmetry. 
Therefore, we demand a zero granular heat flux through 
the circumference r = R which, for T > 0, leads to 



dT 

— (r = R,6,t)=0. 



(8) 



Furthermore, we demand that particle collisions with the 
wall r = R be elastic, and that the wall be impenetrable. 
These lead to slip boundary condition for vg and zero 
boundary condition for v r : 



v r {R,e,t) = 0. 



(9) 



Although somewhat restrictive, boundary conditions (8) 
and (9) allow a lot of insight into the problem, as will 
become clear shortly. The case of no-slip boundary con- 
dition, vg{R,6,t) — 0, is briefly discussed in section III 
D. 

For the narrow-sector-shaped container one should also 
prescribe boundary conditions on the straight segments 
of the sector. However, in the absence of azimuthal vari- 
ations of the hydrodynamic fields, and for vg = 0, there 
is no need of prescribing these boundary conditions ex- 
plicitly. 



B. Global conservation laws 

Equations (2)-(5), together with the "conservative" 
boundary conditions (8)- (9), give rise to two conserva- 
tion laws and a balance relation: 



d_ 

dt 



(P) 
(prvg) 

-pv 2 + pT 



const, 
const. 



a(p 2 t^ 2 ), 



(10) 

(11) 

(12) 



where (...) denotes averaging over the area of the circular 
container: 



(13) 



Equation (10) describes conservation of the total mass 
of the gas. Equation (11) describes conservation of the 
angular momentum of the gas around the center of the 
container. This conservation law is special for the circular 
geometry; it is expected, as neither of the forces in the 
right hand side of Eq. (3) creates torque. Indeed, a 
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pressure gradient cannot create torque, whereas a viscous 
stress causes the momentum only to diffuse but not to 
dissipate. In addition, the wall r = R exerts, according to 
the boundary conditions (9), an effective force on the gas 
only in the radial direction which again does not produce 
a torque. 

Equation (12) describes the global energy loss of the 
gas. As the wall is adiabatic, energy can only escape 
to the internal degrees of freedom of the grains. Equa- 
tion (12) states that, as long as the gas is not "frozen" 
(T 7^ 0), the energy will continue to dissipate. Re- 
lations (10)-(12) are also valid in the sector geometry. 
There Eq. (11) is satisfied trivially, as vg = 0. 



where to — const. For sufficiently small circular con- 
tainers, see below, and if the total angular momentum 
of the gas is zero, the HCS is a global attractor of the 
flow in this geometry. On a qualitative level, the HCS 
stability in small containers is evident from the hydro- 
dynamic equations. If the system size is R, the heat- 
diffusion term in Eq. (5) scales as R~ 2 while the cooling 
term is i?-independent. Therefore, as R becomes small, 
the heat diffusion generates a flow which erases any in- 
homogeneities. 



B. Small-amplitude modes are decoupled 



III. MARGINAL STABILITY OF THE 
HOMOGENEOUS COOLING STATE 

Apart from inertia of the flow, there are two compet- 
ing mechanisms which ultimately determine the spatial 
structure of a freely cooling granular gas. These are 
the inelastic energy loss which drives the clustering- and 
shear-mode instabilities [10, 35], and the heat diffusion 
and viscosity which tend to erase inhomogcneities in the 
gas. The rest of parameters being fixed, the heat diffu- 
sion and viscosity prevail in sufficiently small containers; 
here the granular gas approaches the HCS. When the 
container size is increased, larger-scale perturbations be- 
come possible, and the HCS eventually loses its stability. 
In this section we will find the critical radii of the cir- 
cular container for the clustering modes and the shear 
mode which can bifurcate from the HCS. When the HCS 
is no longer stable, it is the further development of these 
modes which determine the late-time behavior of the gas. 

Linear stability analysis of the HCS was performed by 
McNamara [35], Deltour and Barrat [13], Brey et al. [38] 
and other workers. Here we extend the previous analy- 
sis, that was done in rectangular geometries, to the cir- 
cular geometry. As we are only interested in the bifur- 
cation points (where the system's behavior qualitatively 
changes) and do not need to know the growth rates of lin- 
ear perturbations, we will limit ourselves to a marginal 
stability analysis. 



A. The HCS is stable in small systems 

The HCS is the simplest exact solution of Eqs. (2)-(9). 
Here p = po = const, and v = 0. The energy equation 
becomes 



dT(t)/dt = -A Po T 3 / 2 (t). 



(14) 



and so the gas temperature decays according to Haff's 
law [34]: 



T H (t) 



[Apo(t + t )} 2 ' 



(15) 



Although the HCS is not a steady state of the sys- 
tem in the usual sense, it does become a steady state 
in rescaled variables. The rescaled temperature is 
r(r, t) = T(r, i)/Tf/(t); the density remains unchanged. 
Marginal stability for the gas velocity implies a con- 
stant Mach number: the velocity decays neither faster 
nor slower than the (time-dependent) speed of sound 
y/lTji{t). Therefore, we introduce the rescaled velocity 
u by putting 



v(r,t) = VVM*)u(r,t). 



(16) 



We emphasize that we will only deal with zero- angular- 
momentum flows. Indeed, as the HCS has a zero angular 
momentum, no finite-angular-momentum state can bifur- 
cate from it. 

Let us linearize the hydrodynamic equations (2)- (5) 
around the HCS. We introduce small perturbations to 
each of the variables, 



p = p (l + 5p) 
T = T H (t)(l + 6r) 
V = y/jT H {t)8u, 



(17) 



where p is the average density of the gas. For marginally 
stable states 



d x d , d x n 



(18) 



As a result, the continuity equation (2) becomes 

V • 5u = 0. (19) 
Next, we linearize the energy equation (5): 

2Sp + St - 1 2 k V 2 5t = 0, (20) 



where l K = ^/2k /(Apq) is the critical length for cluster- 
ing instability in a rectangular geometry [10, 35, 36, 38]. 
Finally, we linearize the momentum equation (3): 



Su + l 2 V 2 8u = 



V(S P + St), 



(21) 



where l v = \J 2vq / (kpfy is the critical length for the shear 
instability in a rectangular geometry [10, 35, 38, 39]. By 



■5 



virtue of Eq. (19), the linearized flow is incompressible. 
As a result, the velocity field is purely solenoidal. Now 
we see that the left hand side of Eq. (21) is solenoidal, 
while the right hand side is potential. By the Hclmholtz 
decomposition theorem for vector fields, see e.g. [40], 
each side of the equation must be equal to zero. In other 
words, the shear modes and the density modes (the clus- 
tering modes) decouple. As a result, 

Sp + St = const. (22) 

Now, J D Spd 2 r = by virtue of mass conservation, and 
so JpSrcPr = by virtue of Eq. (20), where D repre- 
sents the circular container. Therefore, Sp = —St: at 
marginality the small density and temperature pertur- 
bations are isobaric. For the velocity we can introduce 
the Stokes' stream function, as the flow is incompress- 
ible. Putting 5u x = dip/dy and 5u y — —dip/dx, we re- 
duce the marginal stability equations to two independent 
Hclmholtz equations 

Sp + l 2 K V 2 Sp = 0, (23) 
V^ + ^V 2 ^ = 0, (24) 

for the clustering modes and shear modes, respectively. 

C. Clustering modes 

The regular independent solutions of Eq. (23) are, up 
to a constant pre-factor, 




FIG. 2: The density field p{r,6) = 1 - 0.1J (r/l K ) for the 
lowest marginally stable clustering mode m = with R = 
Ro = 3.83171... l K . 




cos(m9) 
sm(m9) 



m = 0,1,2,. 



(25) 



where J m (. ..) is the Bessel function of the first kind. 
To satisfy the boundary condition (8), we demand 
J' m {R/l K ) — and obtain a doubly infinite series of solu- 
tions with m = 0, 1, . . . and quantized values of R. The 
first two marginally stable modes are the lowest to = 
mode with R = Rq = 3.83171 ...l K and the lowest m = 1 
mode with R = R 1 = 1.84118 ...l K (see e.g. Ref. [43] 
for the roots of the Bessel functions and their deriva- 
tives). The m = modes are azimuthally symmetric. 
The gas density, corresponding to the lowest of them, is 
p = 1 + CoJo(r/l K ), where the sign of the constant Co 
determines whether the gas is denser at the center or at 
the circumference. An example of this density profile is 
depicted in Fig. 2. For the to = 1 modes the azimuthal 
symmetry is broken; the density profile of the lowest of 
them is p = 1 + C\Ji(r/l K ) cos(#), see Fig. 3. 

As the container radius increases from zero, the first 
bifurcation occurs at R = Ri for the lowest to = 1 mode. 
Therefore, the lowest to = 1 mode, rather than the lowest 
to = mode, is the most unstable clustering mode in this 
geometry. In the following we will sometimes omit the 
word "lowest" in this context. 



D. Shear modes 

Now we proceed to Eq. (24) . The independent regular 
solutions for the stream function are 

j ^){Z(3}'-= ' 1 ' 2 '---- ( 2g ) 

Correspondingly, the velocity components are 



x v -ie±_ mJ ( n ] / -sin(m0) 
= -& = "iWr/M { 



(27) 



sin(m#) 

The boundary conditions (9) for the velocity become 
J m {R/lu) = and J' m {R/l v ) = for to = 1, 2, ... , 

MR/Q =0 for to = 0. 

(28) 

As the first two conditions, for m = 1,2,..., cannot be 
satisfied simultaneously, the only possible type of shear 
modes are the azimuthally symmetric modes to = 0, for 
which the velocity field is 



Su r = and Sug = J\{r/l v ). 



(29) 
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The critical radius for the bifurcation of the lowest m = 
shear mode is R s — 5.13562 ... l v . For a comparison with 
the clustering modes we need to express the critical ra- 
dius for the to = shear mode in units of l K . For a 
gas of hard disks R s = 2.56781 . . . l K , whereas for a gas 
of hard spheres R s = 3.75053 . . . l K . Therefore, one al- 
ways has R s > Ri, and the clustering mode bifurcates, in 
the circular geometry, at a smaller container radius than 
the shear mode. This is in contrast to what happens in 
rectangular geometries, where a shear mode goes unsta- 
ble first [35, 38, 39]. An interesting feature of the shear 
modes in a circular container is that they exhibit zonal 
flows. This is because the flow must have a zero angular 
momentum, as imposed by Eqs. (11), (16) and (18). One 
example is shown in Fig. 4 which depicts the velocity 
field of the lowest to = mode. 



/ s . 
/ / , , 
/ / / ■ , 



N. S. 



-2 



-4 



■ I I I 



-4-2 2 4 

FIG. 4: The velocity field Su = Ji(r/l v )8 of the lowest 
marginally stable shear mode. R — 5.13562... l v . Evident 
is a zonal flow. 

For the no-slip condition at r = R the analysis is 
very similar. Here the critical radius R s for the shear 
mode is determined from equation Ji(R/l v ) = 0. As 
a result, R s = 1. 91585... l K for hard disks, and R s = 
2.79828 . . .l K for hard spheres. Both values exceed the 
critical radius R\ = 1. 84118... l K for inhomogeneous 
clustering. Therefore, for the no-slip boundary, the first 
bifurcation to occur as the system size is increased is 
again the clustering mode bifurcation. Notably, a zonal 
flow is absent for the no-slip boundary. 



IV. NONLINEAR CLUSTERING IN A 
CIRCULAR CONTAINER 

As we have just found, the clustering mode to = 1 is 
the first mode to become unstable as the container radius 
R exceeds R\ = 1.84118 . . . l K . One possible scenario of 
the evolution of the unstable to=1 mode is the forma- 
tion of an inhomogeneous density distribution of the gas 
with no mean flow, v = 0. This state can continuously 
bifurcate from the HCS. In view of Eq. (2) this density 
distribution must be time- independent. In addition, from 



Eq. (3) the pressure gradient must be zero, Vp = 0. The 
density distribution p(r) and the evolution in time of the 
pressure p(t) are governed by the energy equation (6) 
which becomes 



p = -7pV-v+K p 3 ^ 2 V- 



* v(I 

VP VP 



-Ap^y/ 2 , (30) 



with Eq. (8) as the boundary condition. It is convenient 
to rescale the variables: the pressure by its characteristic 
initial value po, the density by the average density po, the 
coordinates x, y by the critical length l K which appears 
in Eq. (20), and the time by the characteristic cooling 
time 



tr. 



The rescaled problem is the following: 



4n = - 2 ^ 1/2 + v • 




dp 
dr 



(r = R) 



0. 



(31) 

(32) 
(33) 



where we abuse notation by denoting the rescaled quanti- 
ties by the same symbols as the physical ones. Averaging 
Eq. (32) over the (rescaled) area of the container, we ob- 
tain a global balance equation for the pressure 



P 



p 



,3/2 



= -2<p 1/2 > 



(34) 



where (...) denotes averaging over the area of the circu- 
lar container. Thus, the time-independent density pro- 
files are determined from the nonlocal problem 



1 

3 

dp 
dr 



vV 3/2 + (p 1/2 ) - p 1/2 = o, 

(r = R) = 0. 



(35) 



This problem is a generalization to 2D of the similar 
problem for a narrow-channel flow [36]. A 3D version 
of this equation was derived by Fouxon [37]. He dis- 
cussed spherically-symmetric solutions with a single den- 
sity maximum at r — R. In 2D these correspond to the 
lowest to = mode, see Fig. 2. As the to = 1 mode 
actually goes unstable first in a circular container, any 
azimuthally-symmctric solution of Eq. (35) must be un- 
stable with respect to perturbations breaking this sym- 
metry. We checked that a similar feature persists in 3D: 
the first mode to become unstable as the radius of a 
spherical container is increased from zero is not spher- 
ically symmetric. As a result, in a spherical container, 
the spherically-symmetric solutions reported in Ref. [37] 
must be unstable with respect to perturbations break- 
ing this symmetry. The azimuthally- and spherically- 
symmetric solutions (in 2D and 3D, respectively) be 
come, however, relevant in narrow-sector (correspond- 
ingly, narrow-cone) shaped containers, where heat dif- 
fusion erases any non-radial structure. We will deal with 
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the narrow-sector geometry later on, in section V. In 
this section we continue to deal with a circular container. 
First of all, we will perform a weakly nonlinear analysis 
of Eq. (35) to determine the character of bifurcation of 
the HCS. As we will see, the bifurcation at hand is a 
pitchfork bifurcation, where the state parameter (a den- 
sity contrast defined below) grows as a square root of 
the difference between the container radius and the crit- 
ical radius. A pitchfork bifurcation can be either super- 
critical, or sub-critical [42], thus determining the stabil- 
ity of states below and above the bifurcation. Since this 
important characteristic is determined by the sign of a 
numerical coefficient, an accurate calculation is needed 
to ascertain it. 



A. Bifurcation analysis 

Equation (35) can be simplified by a transformation of 
the unknown function p and of the rescaled radial coor- 
dinate r: 



A~ 2 r 



(36) 



where A = (p 1 ^ 2 )- We obtain a parameter-free nonlinear 
Poisson equation with a Neumann boundary condition at 
f = R= A 2 R: 



-V 2 $ + 1 - $~ 1/3 



0, 



(37) 



Once $ is found, A can be determined from A 2 = 
($~ 2 / 3 )f, where the averaging is over the circular con- 
tainer of rescaled radius R. The latter relation follows 
from the condition (p) = 1 for the rescaled density p. 

In the new coordinate f , the critical radius for insta- 
bility of the mode m = 1, found in section IIIC, is 
Ri = R\ = 1.841 Let us consider a circular con- 
tainer with a radius slightly above or below R\ and seek 
for a weakly- nonlinear solution of Eq. (37). At the risk of 
irritating the reader, we will transform the radial coor- 
dinate one more time, by introducing \ = {Ri/R)?-> and 
also define the parameter k = R/ R\ which will serve as 
the eigenvalue of the nonlinear eigenvalue problem 



iv 2 $ + £: 2 (l-$~ 1/3 ) 







j-(Ri,e) = o. 



(38) 



Now we expand both the eigenfunction and the eigen- 
value in a series employing a perturbation amplitude e as 
a small parameter: 



$(X,0) = l + e Vl (x,9) + e 2 ^ 2 ( X ,e) 
+ e\ 3 ( X ,6) + ... 



l + ek 1 + e 2 k 2 



(39) 



where, according to our marginal stability analysis, ipi — 
J\(x) cos(0). In the second order in e we obtain 

\7 2 tp 2 +<f2 = - 2fci<pi 



(40) 



The solvability condition of the inhomogeneous linear 
equation is the orthogonality of its right hand side to the 
function <fi\. This yields fci = 0. Solving the resulting 
equation by a Fourier decomposition, we find 

^2-^Co(x) + Jc 2 (x)cos(20), (41) 

where the radial functions (o and £2 are solutions of the 
forced Bessel equations 



0,2, 



(Km 

dx 



0. 



whereas 



1 d d 

xdx X dx 



1 



(42) 



(43) 



is the Bessel operator of order m. 
In the third order in e we obtain 

^2 4 14 3 „, 



0. 



(44) 



The solvability condition again demands to choose k 2 so 
that the right-hand side be orthogonal to ipi. This is 
achieved by setting 



k 2 = 



-0.02153.... 

(45) 

Here the integration J J D ■ ■ ■ xdxdd is carried out over 
the circular container D of radius Ri . With this value of 
k 2 we can calculate the dependence of e on SR — R — R\ : 

''' " 1 1 -SR. (46) 



k 2 



Rik 2 



As k 2 < 0, SR can only take negative values. That is, 
somewhat surprisingly, a time-independent solution de- 
scribing an inhomogeneous density distribution is pos- 
sible only in a container with radius smaller than R\. 
Going back to the rescaled variables (that is, two trans- 
formations back), we find 



1 + 



(Mx) 2 ) 



18 



x 2 

C 1 



SR = R — R\ 



I8k 2 



SR, 



p{r,9) - 1 \eJ, (^r) ■ 



(47) 
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Now we can determine the form of the bifurcation di- 
agram in the (rescaled) original coordinates. Taking 
5R = R — R\ to be the bifurcation parameter, and the 
density maximum p max to be the bifurcation amplitude, 
we obtain 



Pmax =s 1 + 3.1432.. .y/R 1 - R, R x - R < i2 x . (48) 

This is a subcritical bifurcation, where inhomogeneous 
solutions exist only at R < R\. The time-independent 
inhomogeneous solutions form a continuous set which can 
be obtained by rotating one of the solutions by an arbi- 
trary angle around the center of the container. Figure 5 
presents the bifurcation diagram alongside with its exten- 
sion obtained by numerical integration of Eq. (38) (using 
the pdenonlin algorithm of the Matlab PDE toolbox) 
and transformation back to the (rescaled) original vari- 
ables. The well-known properties of supercritical bifurca- 
tions in ID [42] make it possible to determine the stability 
of the branches shown in Fig. 5 without further calcula- 
tions. The arrows in Fig. 5 schematically show the di- 
rection of the evolution of the flow depending on whether 
the branch is stable or unstable. We found that the max- 
imum density of the unstable states diverges at a finite 
subcritical container radius, R ~ 1.7795 ~ 0.9666 R\, 
see Fig. 5b. This suggests that this is a critical radius 
separating, in a proper space of parameters, the flows to 
those which approach the HCS from those which do not. 
Figure 6 shows the numerically calculated density map 
of the "border state" corresponding to this critical radius 
of the container. 

One salient feature of the bifurcation diagram in Fig. 
5 is that there are no stable time-independent inhomoge- 
neous solutions at R > R\. This result is in a strik- 
ing contrast with that for the narrow-channel geome- 
try, where the pitchfork bifurcation for the clustering 
mode is super-critical, and where the energy loss and 
heat diffusion balance each other to produce stable time- 
independent inhomogeneous states [36] . What is then the 
ultimate state of a freely cooling granular gas at R > R{! 



B. Hydrodynamic simulations in a circular 
container 



Our first step in answering this question was perform- 
ing a series of 2D hydrodynamic simulations in this geom- 
etry, see Appendix A for the computational details. Af- 
ter a proper rescaling the time-dependent hydrodynamic 
equations include only one dimensionless parameter, pro- 



portional to 1 



see Ref. [36] and Appendix A. We 



chose the restitution coefficient a 
corresponds to e\ — k A/2 = 1 — 



0.94868... which 
2 = 0.1. The hy- 
drodynamic equations, subject to boundary conditions 
(8) and (9), were solved with the (zero- viscosity) Vul- 
can/20 code, see Appendix A. An Eulcrian numerical 
scheme on a non-uniform polar mesh was used. Here 
we present the results of a typical simulation, performed 
in (one half of a) circular container of rescaled radius 



1.3 = 



1.2 



1.1 



X 

I \ t 



1.836 1.838 



R 



1.84 1.842 




FIG. 5: (Color online.) Time-independent inhomogeneous 
cooling states in a circular container exhibit a sub-critical 
pitchfork bifurcation. Shown is the gas density at the cir- 
cumference r = R versus the container radius R. The solid 
lines denote stable states, the dotted lines denote unstable 
states. The horizontal line p m ax = 1 depicts the HCS. The 
dashed lines depict analytical solutions. The circles show re- 
sults of numerical integration of Eq. (38). All the quantities 
are plotted in rescaled units. The arrows schematically show 
the direction of the evolution of the flow depending on the 
stability properties of the branches. Panel (b) is an extension 
of panel (a) to smaller container radii. 



R = 2.3. This value is larger than the rescaled critical 
radius Ri = 1.84 . . . for the m = 1 clustering mode, but 
smaller than the critical radii for the shear mode and for 
all other clustering modes in this geometry. The initial 
conditions were 



p(r,9,t = 0) = 1 + J 1 (R 1 r/R)cos{9), 
v(r,6,t = 0) = 0, 
p(r,9,t = 0) = 1. 



(49) 



The simulation results are illustrated in Figs. 7 and 
8. Figure 7 shows density map snapshots at different 
(rescaled) times. One can see that the gas develops a 
density peak at the point r = R, 6 = Q. Figure 8 shows 
the evolution of this density peak in time. The numerical 
run was stopped when the (rescaled) maximum density 
reached about 300, and a considerable amount of gas 
was already concentrated in a region comparable with 
the mesh size. One can see that, until this time, the 
density growth accelerates in time. These results suggest 
that the gas is developing thermal collapse: a finite-time 
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FIG. 6: The density map of the unstable inhomogeneous 
steady-density state in a circular container of radius R ~ 
1.7795 ~ 0.9666 Ri. This state, obtained by solving numer- 
ically Eq. (38), serves as a "border" between those cooling 
flows which approach the HCS and those which do not. 
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FIG. 7: The density history of a freely cooling granular gas in 
a circular container of rescaled radius R = 1.3-Ri as obtained 
in a hydrodynamic simulation for initial conditions (49). One 
half of the container is shown. The (rescaled) times of the 
snapshots are t = 1.7 (a), 19.1 (b), 32.12 (c) and 39.45 (d). 



density blowup which, contrary to what happens in a 
narrow-channel flow [36], is not arrested by heat diffu- 
sion. To substantiate this hypothesis, we will now switch 
to a narrow-sector geometry, see Fig. 1. 
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FIG. 8: The maximum gas density versus time corresponding 
to Fig. 7. 



V. TIME-INDEPENDENT INHOMOGENEOUS 
COOLING STATES IN A NARROW SECTOR 

In a narrow-sector-shaped container the heat diffusion 
suppresses any flow structure in the azimuthal direction. 
This setting, therefore, defines a purely radial (ID) flow 
while retaining an essential 2D feature: polar metrics. 

As the (azimuthally-asymmetric) m = 1 clustering 
modes are now suppressed, the first mode that goes un- 
stable, as the sector radius R is increased from zero, is 
the lowest m = clustering mode, see Sec. IIC. Let the 
sector radius be slightly larger than the critical radius 
i?o = 3.83171Z K for the lowest m — mode. The unsta- 
ble m — mode will start to grow. What is the long- 
time behavior of the gas in this geometry? To answer 
this question, we will first identify a time-independent 
inhomogeneous solution: a possible candidate for the de- 
scription of a long-time behavior. We will show that the 
bifurcation of the HCS is trans-critical here, and deter- 
mine its sign by performing a weakly nonlinear analysis. 



A. Bifurcation analysis 

Let the rescaled sector radius be slightly larger, or 
slightly smaller, than Rq = 3.83171 .... We again use the 
transformed equation (37). Dropping the ^-derivatives, 
we obtain the nonlinear eigenvalue problem 



l rf 2 $ 
3 rff 5 " 



1 d<$> 

3f df 



- 1 - $ 



-1/3 



0. 



R)=0, 



(50) 



the sector radius serving as the eigenvalue. As the HCS 
corresponds to $ = 1, we substitute <f> = 1 + eJ (r) + 
e 2 (p2(r) in Eq. (50) and neglect terms of orders e 3 and 
higher. The resulting equation is a forced Bessel equa- 
tion: 



</?2 



df 2 



1 dip 2 
f df 



(51) 
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This equation can be solved via variation of parameters, 
and we obtain 



$ = 1 + eJo(r) - -eVo(f) / J (f'YY (r')f'df' 

Jo 

+| e 2 r (f)^ r Mf'ff'd? , (52) 

where Yq(. . .) is the Bessel function of the second kind. 
This inhomogeneous solution must hold at R ^ Rq, so 
e depends on R. We can seek R = Rq + 5R, where 
SR <C Rq. Forcing the solution to obey the boundary 
condition 



— (Ro + SR) = 0, 
ar 



(53) 



we obtain, in the lowest order, a linear relation between 
e and 5R: 



5R = 2.223... SR. (54) 



Returning to the rescaled variables r and p, we obtain 
the density at the arc of the sector, r = R: 



p{r = R) = l--J (Ro)e = 1+0.5968... (R-Ro), (55) 



to first order in R — R$. 
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FIG. 9: Inhomogeneous cooling states below and above the 
bifurcation point R = Rq, obtained by numerical integration 
of Eq. (50) for different sector radii. The solid curves show 
stable solutions, the dashed curves show unstable ones. The 
quantities are plotted in rescaled units. 

It is clear from this calculation that, at least for system 
sizes close to Rq, the cooling gas has two possible steady- 
density states. The first is the HCS p — 1. The other is 
an inhomogeneous state where, for R < Rq, the density 
is peaked at the vertex of the sector, and for R > Ro 
the density is peaked at the arc. Several examples of 
density profiles of these inhomogeneous states, obtained 
by numerical integration of Eq. (50), are shown in Fig. 
9. 




FIG. 10: (Color online.) Time-independent inhomogeneous 
cooling states in a narrow-sector-shaped container exhibit a 
trans-critical bifurcation. Shown is the gas density at the arc 
r = R versus the sector radius R. The solid lines denote sta- 
ble states, the dashed lines denote unstable states. The HCS 
is represented by the horizontal line p(R) = 1. The inclined 
straight line is the relation (55). The results of numerical in- 
tegration of Eq. (50) are shown as circles. The arrows indicate 
the conjectured flow in a system's phase space. All the quan- 
tities are plotted in rescaled units. Panel (b) is an extension 
of panel (a) to larger and smaller sector radii. 



Equation (55) describes a trans-critical bifurcation [42] 
in which two states "collide" and exchange their stabil- 
ity. The bifurcation parameter is SR. As R increases, 
the HCS remains stable until R = Rq. Here the HCS 
"surrenders" its stability to the inhomogeneous state. 
Correspondingly, the inhomogeneous state is unstable at 
R < Ro- The bifurcation diagram of the system is shown 
in Fig. 10 a. 

The bifurcation diagram can be extended beyond 
the weakly-nonlinear regime by numerical integration of 
Eq. (50). Figure 10 b shows the computed stable branch 
of the inhomogeneous solutions at R > Rq, and the un- 
stable branch at R < Rq. Note the exponential depen- 
dence, at large radii, of the maximum density on the 
sector radius. It has the same origin as the exponential 
system-size dependence of the maximum density of the 
inhomogeneous cooling state in a long narrow channel 
[36]. Our numerics show that the unstable branch ends 
at R ~ 3.25, where p(R) — 0.77. That is, there are no 
time-independent inhomogeneous solutions for the den- 
sity at smaller radii. Figure 12 shows, for the unstable 
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p{R) = Cpi{R;R) + 1 - C, so 



FIG. 11: A steady-state density profile of a freely cooling 
granular gas in a narrow-sector-shaped container of rescaled 
radius R — 7.57. The quantities are plotted in rescaled units. 



C 



Pi{R\R)- 
p{R) - 1 



(57) 



For initial conditions located above the unstable branches 
of Fig. 10, the gas will approach the corresponding sta- 
ble state: the HCS at R < Rq, and the inhomogeneous 
cooling state at R > Rq. There are no stable states, how- 
ever, when starting below the unstable branches, see Fig. 
10. This suggests that, as in the circular container, the 
cooling gas develops thermal collapse. In the following 
section we will present analytical and numerical evidence 
supporting this conjecture. 



VI. THERMAL COLLAPSE IN 2D 
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FIG. 12: The gas density at the vertex of a narrow sector 
versus the sector radius for unstable inhomogeneous solutions 
at R < Ro- The quantities are plotted in rescaled units. 



branch, how the density at the vertex, p(r — 0), grows 
when R goes down. 

When increasing the sector radius, the density profile 
becomes more and more peaked at r = R, so one can call 
these states localized [37]. In Fig. 11 one can see the 
radial density profile for R = 7.57. 

The bifurcation diagram in Fig. 10 suggests a "phase 
diagram" of the cooling gas in the sector geometry. Let 
us denote the time-independent inhomogeneous solution 
of Eq. (35) in a sector of radius R by pi{r; R). Consider 
a continuous one-parameter family of initial conditions 
such that the density profile is a linear combination of 
the inhomogeneous solution pi(r \ R) and the HCS p = 1: 



p(r,t = 0) = Cpi{r;R) + l 
v(r,t = 0) = 0, 
p(r,t = 0) = 1. 



(56) 



The constants are chosen so that (p) — 1 . Let us identify 
an arbitrary point [R, p] of the plane of Fig. 10 with this 
initial condition. Then the constant C is determined by 



A. Reduced ideal hydrodynamics close to 
singularity 

That a freely cooling dilute granular gas exhibits ther- 
mal collapse, in any dimension, within the framework of 
ideal granular hydrodynamics, is by now well established 
[24, 29, 30, 36, 37]. In this subsection we will derive a set 
of reduced ideal hydrodynamic equations which describe 
thermal collapse in 2D. In the next subsection we will 
obtain scaling laws, in time and in space, which charac- 
terize the thermal collapse. Then we will check whether 
heat diffusion and other processes, that we will neglect, 
can become relevant near the singularity. 

Assuming that thermal collapse does happen, and that 
heat diffusion and viscosity are irrelevant in the vicinity 
(both in time, and in space) of the density singularity, one 
can considerably simplify the hydrodynamic equations 
Eqs. (2), (3) and (6). In the light of results of Refs. 
[24, 29, 30, 36, 37], an important additional simplification 
comes from the assumption that the gas pressure in the 
collapse region is approximately uniform. In a formal 
language, see e.g. Ref. [41], we eliminate the acoustic 
modes by postulating a perturbation expansion for the 
gas pressure, 



(58) 



limit ourselves to the leading-order approximation, p — 
p(°'(t), and discard the momentum equation (3). Finally, 
we neglect both terms in the left hand side of Eq. (6), as 
they remain bounded at the density singularity, and keep 
only terms which blow up at the singularity. We arrive 
at the following reduced equations: 



V ■ (pv) = 0, 



dp 

i + 

7P7 1/2 V • v 



-hp 



1/2 



(59) 
(60) 



where p s is the gas pressure at singularity. Equation (60) 
describes a balance between the adiabatic compression 
heating and the inelastic cooling. 
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Set of equations (59) and (60) is not closed, as there 
are only two scalar equations for the density and the two 
velocity components. To obtain an additional equation, 
for the velocity, one should go to the subleading order in 
Eq. (58), as it was done in Rcf. [41] in the case of a 2D 
gas flow driven by heat diffusion. Instead, we will make 
another simplifying assumption: that the gas flow near 
thermal collapse is (at least locally) azimuthally sym- 
metric. Note that, if thermal collapse develops at the 
vertex of a narrow-sector-shaped container, this assump- 
tion is satisfied trivially. For an azimuthally-symmetric 
flow Eqs. (59) and (60) make a closed set. 



B. Scaling laws of thermal collapse 

Let us rescale t by t s = 2^/(Kp 1 J 2 p^ 2 ), p by the av- 
erage gas density p , r by l K and v by l K /t s . Then, for a 
<i-dimcnsional thermal collapse (azimuthally symmetric 
in 2D or spherically symmetric in 3D) Eq. (60) acquires 
a parameter- free form: 



d 



Or 



(r d -\)=-2p^. 



(61) 



where r is a local radial coordinate with the origin at the 
point of developing singularity. Substituting Eq. (61) in 
Eq. (59) we obtain 



(62) 



This equation becomes very simple in Lagrangian coor- 
dinates which transform the full time derivative on the 
left-hand-side into a simple time derivative. We can de- 
fine a Lagrangian mass coordinate 



m(r, t) 



f ' p(r',t)r' d - 1 dr' 
Jo 



(63) 



being the metric in d dimensions. Up to a d- 
dependent factor, m(r, t) is the mass of the gas contained 
in a d-dimcnsional sphere of radius r around the point of 
developing singularity. After the transformation Eq. (62) 
becomes remarkably simple: 



dw i 



= -1, 



Qj- \m— const 

where w(m,t) = p^ 1 / 2 {m,t). The solution is 
w{m, t) = —t + wo(m), 



(64) 



(65) 



where wo (to) is the initial condition (with a minimum at 
to = 0) which should be chosen at time when the reduced 
set of equations (59) and (60) already holds. The time of 
singularity is therefore t* = u> (0). The scaling laws near 
collapse are determined by the behavior of the function 
Wo(m) at small m. To elucidate this behavior, let us re- 
turn to the density p in the original Eulerian coordinates. 



Here, a generic density profile near the singularity has a 
quadratic dependence on r: 



p (r) ~ a-br 2 



(66) 



where a > and b > 0. Using this local profile we can 
calculate, for this "initial condition" , the local asymptote 
m(r) near r = 0: 



m(r) 



po(r')r' dr' 



in the leading order in r. Therefore, 

2/d 



po(m) ~ a - b ( - 



to 



2/d 



This follows 



w {m)~a- l/2 + cm 2/d , 



(67) 



(68) 



(69) 



where c > is a constant. Now we see that the den- 
sity blowup occurs at t* = a -1 / 2 . Now we return to 
Eq. (65) and observe that, sufficiently close to the singu- 
larity point, the solution is self-similar: 



w 



(to, t) ~ U - t + cm 2/d = t W{m/T d/2 ) , (70) 



where t — t* — t and W(z) = 1 + cz 2 l d . Transforming 
back to the Eulerian coordinate, we obtain 

d A r dm ' a r a * 2 w2 ( m' \ 

r =d L — , 'L dmT2w \— 2 ) 

= - 2+d/2 ^)' ^ 

where J r (. . . ) is a shape function. Equations (70) and 
(71) yield the dynamical length scale of the ideal singu- 
larity 



*~(t„-t)", 



1/2 + 2/d. 



(72) 



For d = 1 (the narrow-channel geometry) we obtain (3 — 
5/2, in agreement with the result of Fouxon et al. [29, 30], 
whereas for d = 2 and 3 we obtain new exponents /3 = 3/2 
and 7/6, respectively. The density asymptote close to the 
singularity can be written as 



P(r,t) 



(73) 



where £ = r/r 13 is the similarity variable, and ^ is a 
shape function such that \P(0) = 1. The density blows 
up at r = as p = (t r — t)~ 2 independently of dimension 
[37]. 

At large distances, r » £ (but still much smaller than 
an external length scale of the flow) the density exhibits 
a time-independent power-law tail. Indeed, at m > T d l 2 
asymptote (70) becomes w(m,t) ~ cm 2 ^ d . Then, using 
Eq. (71), we obtain r d ~ m 1+i / d . As a result, 



(74) 



13 



and 



p(r » t) <~ r v , i] = — 



Ad 
4 + rf 



2 



(75) 



For d = 1 this yields = —4/5, in agreement with Refs. 
[29, 30], whereas for d = 2 and 3 one obtains new expo- 
nents rj = —4/3 and —12/7, respectively. 

To find the scaling relations for the gas velocity we 
integrate Eq. (61) 



' Jo 



(76) 



and change to the self-similar coordinate £: 



"V(0- 



(77) 

When £ ^> 1, the velocity field exhibits a time- 
independent power-law tail, V <~ — where /x = 1 — 1/(3. 
Again, for d = 1 we recover the previous result p, = 3/5 
[29, 30]. In 2D /z= 1/3, in 3D // = 1/7. 



C. Can heat diffusion arrest thermal collapse? 

Does heat diffusion, unaccounted for in ideal descrip- 
tion, arrest the collapse? In ID the answer to this ques- 
tion is affirmative [36]. What happens in higher dimen- 
sions? Using the scaling relations in the vicinity of ther- 
mal collapse, that we have just found, we can estimate 
the ratio of the heat diffusion term to the cooling term. 
In the original units of Eqs. (2)- (6) we obtain: 



«oV • Wp/p V(p/p)] 



JLL 
Po L 



t. 



(78) 



One can see that the case of d = 1 is dramatically differ- 
ent from the cases of d = 2,3. For d = 1 heat diffusion 
cannot be neglected when time is sufficiently close to the 
t*. As the heat diffusion is a stabilizing factor, one can 
argue that, in this case, it may arrest thermal collapse. 
Indeed, it was shown in Ref. [36] that, in the narrow- 
channel geometry, heat diffusion ultimately balances the 
unstable density growth, and a time-independent density 
profile sets in. On the contrary, for d = 2 or 3 the heat 
diffusion remains irrelevant and cannot arrest the density 
blowup. 

Estimate (78) implies that, in a narrow-sector-shaped 
container, thermal collapse is only possible exactly at 
the vertex and not in any other location. Indeed, let 
r = r s > be such another location. The differential of 
the local Lagrangian coordinate, measured from r = r s , 
is proportional to dr. This setting is essentially ID, and 
so heat diffusion will arrest the collapse [36]. On the 



contrary, the differential of the Lagrangian coordinate, 
measured from r — 0, is proportional to rdr. That is, it 
is essentially 2D, and the scaling laws of the ideal singu- 
larity persist in spite of the heat diffusion. 

In Appendix B we checked the consistency of other 
assumptions behind the reduced ideal set of equations 
(59) and (60). As it turns out, one assumption - about 
the (almost) uniform pressure profile - breaks down at 
t very close to t s . We will briefly discuss this issue in 
section VII. 



D. Hydrodynamic simulations in a narrow sector 

To test our theoretical predictions for a narrow-sector 
geometry, we solved numerically the rescaled hydrody- 
namic equations, see Appendix A. We assumed that the 
hydrodynamic fields depend only on the radial coordi- 
nate, whereas vg — 0, and used the zero- viscosity Vul- 
can/ ID code, see Appendix A. Here we will describe two 
different simulations for R = 1.3Rq (above the critical 
radius for the lowest m = mode) . The restitution coef- 
ficient a — 0.994987 . . . was chosen to satisfy the relation 
KoA/2 = 1 — a 2 = 0.01, whereas 7 was taken to be 2. The 
first initial condition describes an isobaric configuration 
with a density peak at r — R: 



p(r,t = 0) = 1-0.3 



4 /7rr \ 



v(r, t 
p(r, t 



0) 
0) 



0, 
1. 



(79) 



Here we should expect, see Fig. 10, that the density 
profile will approach the corresponding time- independent 
inhomogeneous solution that we found in section V. This 
is indeed what the simulated dynamics shows, see Fig. 
13. 

The second initial condition differs from the first in 
that the initial density is now peaked at the vertex r = 0: 



p{r,t = 0) = 1 + 0.3 



4 /7rr\ 

^ + C0S {r) 



v(r,t = 0) 
P(r,t = 0) 



0. 
1. 



(80) 



The bifurcation diagram, see Fig. 10, suggests that the 
gas will develop thermal collapse. The simulated history 
of the gas density at the vertex is depicted in Fig. 14. 
Shown is the inverse square root of the density versus 
time. This time dependence is in excellent agreement, 
until the latest simulation times we could reach, with the 
finite-time density blowup p(r = 0,t) <~ (t* — t)~ 2 . Next 
we turn to the spatial behavior of the evolving density 
profiles. Figure 15 shows, in a log-log scale, several den- 
sity profiles at different times. It is evident that, as the 
gas density blows up at the vertex, the density profile de- 
velops a time-independent power-law tail with exponent 
—4/3 as predicted by our reduced ideal theory. 
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FIG. 13: (Color online.) The density profile p(r, t) approaches 
a time- independent inhomogeneous state in a sector of radius 
R — 1.3-Ro- The profiles were obtained in a hydrodynamic 
simulation with initial conditions (79) at times t = (line 
1), 500 (line 2), 2 ■ 10 4 (line 3) and t = 10 9 (line 4). The 
unnumbered line shows the stationary profile obtained from 
the solution of Eq. (50). The density p is measured in units 
of po, and r is measured in units of l K . 
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FIG. 14: (Color online.) Inverse square root of the gas density 
at the vertex of a narrow sector, versus time, obtained in a 
hydrodynamic simulation with initial conditions (80). The 
density blows up as ~ (t* — t)~ 2 . The dashed straight line is 
a guide for the eye. 



A more detailed check of the density blowup addresses 
the coefficient of the scaling relation at r — as pre- 
dicted by the ideal theory. p(0, t) should behave, close to 
blowup, as 



p(0,i)-> 



7 



p s (U - tf 



(81) 



where the time is measured in units of t c , the pressure in 
units of the initial pressure po, and the density in units of 
Pq. Using data from the simulation for the initial condi- 
tions (80), we plotted the quantity 7p(0, i)~ 1/2 p(0, t)~ x / 2 
versus time. This plot is shown, in a log-log scale, in Fig. 
16. One can see that the curve exhibits a power law 
with exponent —2 and passes through the point (1,1), 
showing excellent agreement with Eq. (81) in both the 



FIG. 15: (Color online.) Development of a density blowup 
at the vertex of a narrow sector of radius R = l.SRo in a 
hydrodynamic simulation with initial conditions (80) . Shown 
are the density profiles at times t — (line 1), 117.35 (line 
2), 122.81 (line 3), 123.16 (line 4) and 123.18 (line 5) in units 
of the cooling time t c . The density is plotted in the units 
of the average density po. The time-independent part of the 
observed profiles exhibits a power law with exponent —4/3 
predicted by our theory. 



power law and the coefficient. Now we proceed to ver- 




t, -t 

FIG. 16: (Color online.) Shown is jp(0, t)- 1/2 p(0, t)~ 1/2 ver- 
sus t* — t when approaching the density blowup for the hy- 
drodynamic simulation with initial conditions (80). The time 
of the singularity i* = 123.19 was the only fitting parameter. 
The dashed line shows the function (t* — t)~ 2 . The time is 
measured in units of t c . 

ify the self-similarity of the density profiles at different 
times. Figures 17 depicts the shape function ^(^) in- 
troduced in Eq. (73). Notably, the validity range of the 
self-similar asymptote expands with time. We also veri- 
fied the power-law behavior in time of the dynamic length 
scale i(t). As a measure of £(t) in the simulations we de- 
fined the radial coordinate L(t) where the gas density was 
equal to one half of the density at r = 0. Figure 18 shows 
the time dependence of L(t), and excellent agreement 
with our analytical prediction in 2D, £(t) ~ L(t) ~ t 3 / 2 
is observed. 
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FIG. 17: (Color online.) Emergence of the self-similar shape 
function of the density for the hydrodynamic simulation with 
initial conditions (80). The times are 108.23 (line 1), 120.8 
(line 2), 122.8 (line 3), 123.12 (line 4), 123.18 (line 5) and 
123.18 (line 6), measured in units of t c . p s ~ 7.765 x 10 . 




FIG. 18: (Color online.) The dynamical length scale of ther- 
mal collapse versus the time to singularity for the hydrody- 
namic simulation with initial conditions (80). The dashed line 
is shown to guide the eye. 



Finally, most of our analytical results were obtained 
on the assumption that the gas pressure in the vicin- 
ity of thermal collapse remains uniform. We tested this 
assumption directly by following the simulated pressure 
profiles at different times. Figure 19 shows several pres- 
sure profiles in the whole sector, and a closeup in a vicin- 
ity of r = 0. One can see that, as thermal collapse pro- 
gresses, the pressure remains close to uniform in the ex- 
ternal region but exhibits some steepening near r — 0. 
As the gas density varies by several orders of magnitude 
during these times, the spatial and temporal variations of 
the pressure are insignificant. At the end of the simula- 
tion, the spatial length scale of the pressure is still much 
larger (by more than two orders of magnitude) than the 
dynamic length scale l(t). That is, the pressure remains 
close to uniform in the (shrinking with time) collapse 
region. Still, the pressure steepening may signal a fu- 
ture breakdown of the uniform-pressure approximation 
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FIG. 19: (Color online.) The pressure profiles at t = 121.14 
(line 1), 121.96 (line 2), 122.84 (line 3) and 123.17 (line 4) 
for the hydrodynamic simulation with initial conditions (80). 
The inset shows a close-up of these profiles near r — 0. As the 
blowup time approaches, the pressure profile steepens near 
r = 0. The steepening is very modest, however, compared 
with the giant density growth at these times. Notice also 
that the characteristic length scale of steepening is larger by 
more than two orders of magnitude than the dynamic length 
scale £(t) of thermal collapse. 

at times extremely close to collapse, see Appendix B. 



VII. DISCUSSION 

We have found that, in a circular container, the critical 
radius of the container for the clustering mode instabil- 
ity of a dilute nearly elastic granular gas is smaller than 
that for the shear mode instability. This is in contrast 
with rectangular geometry, where the shear mode bifur- 
cates at a smaller system size than the clustering mode 
[10, 35, 38, 39]. This fact enabled us to investigate a 
strongly nonlinear dynamics of a linearly unstable clus- 
tering mode in a 2D setting, without worrying about the 
shear mode. We have presented analytic and numeric ev- 
idence for thermal collapse, unarrested by heat diffusion, 
in a freely cooling dilute granular gas in 2D. Our results 
imply that Navier-Stokes hydrodynamics of dilute gran- 
ular gases exhibits, in dimension higher than one, finite- 
time density blowups. This is in striking contrast to the 
"wave breaking" singularities of an ideal flow (both in a 
"classical" gas [4], and in a granular gas [29, 30]) which 
are arrested by the viscosity and heat diffusion, bringing 
about smooth shock-wave solutions. 

We have derived the dynamic scaling laws of ther- 
mal collapse by assuming a locally uniform pressure pro- 
file. All the scaling laws, and the uniform-pressure sce- 
nario, have been corroborated in our hydrodynamic sim- 
ulations. An order-of-magnitude estimate, presented in 
Appendix B, suggests that the uniform-pressure scenario 
must break down sufficiently close to the singularity. We 
have not observed this breakdown in our hydrodynamic 
simulations: the uniform-pressure scenario agreed very 
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well with the simulations until the latest simulation times 
we were able to reach. This may be related to some 
small numerical factors which improve the applicability 
of the uniform-pressure scenario beyond what our order- 
of-magnitude estimate predicts. Future work should ad- 
dress this issue and find out whether the uniform-pressure 
scenario gives way to the zero-pressure scenario where a 
different type of singularity - that of a free flow - devel- 
ops [22, 28]. 

It would be interesting to investigate the hydrodynam- 
ics of clustering in a free cooling granular gas in 3D. One 
can speculate that the 2D singularities, analyzed in this 
work, are the generic type of singularity in 3D. Indeed, a 
locally-2D pinching may be more common than a fully- 
3D collapse. 

The dynamic scaling laws, that describe the den- 
sity singularity, provide a sharp characterization of a 
strongly nonlinear stage of clustering instability. When 
the growing-in-time gas density peak reaches a fraction 
of the close-packing density ~ <J~ d , the dilute hydro- 
dynamic equations break down. Our solutions describ- 
ing thermal collapse represent, therefore, an intermedi- 
ate asymptotic of the true flow. Importantly, the du- 
ration of this intermediate asymptotic can be made ar- 
bitrary long if the average density of the gas is made 
sufficiently small compared with the close-packing den- 
sity. Furthermore, the appearance of singularities in the 
dilute-hydrodynamic solutions signals the formation of 
closely-packed clusters of particles. The details of forma- 
tion and dynamics of the closely-packed regions are ob- 
viously beyond a dilute-gas hydrodynamic description. 
A satisfactory description of granular gases at moder- 
ate densities (at a fraction of the close-packing density) 
is provided by a Navier-Stokes granular hydrodynamics 
with the equation of state of Carnahan and Starling [44] 
and semi-empiric transport coefficients derived from ki- 
netic theory in the spirit of Enskog approach [45] . A hy- 
drodynamic description at even higher densities is even 
more challenging, although some particular flow regimes 
at very high densities have been successfully described 
by using empirical constitutive relations [12, 46-48] and, 
when necessary, by taking into account a multi-phase 
character of the flow [49, 50]. 



Appendix A. Hydrodynamic simulations 

We solved numerically a rescaled and transformed ver- 
sion of hydrodynamic equations Eqs. (2)-(5). The rescal- 
ing reduces considerably the number of relevant param- 
eters, whereas the transformation dramatically reduces 
the computation time. 

The rescaling procedure is essentially the same we used 
at the beginning of section IV: the coordinates x, y are 
rescaled by the characteristic length scale l K which ap- 
pears in Eq. (20), the time by the characteristic cool- 
ing time t c = 2/(ApJ /2 pJ /2 ), the gas velocity by l K /t c , 
the density by the average density po, and the temper- 
ature by its typical initial value Tq. The rescaling does 
not change the continuity equation (2). The momentum 
equation (3) becomes 



dv 



+ (v • V)v 



= -V(pT) + e 2 V • (VTS), (Al) 



where t\ = KoA/2, and e 2 = ^oA/2. As we assume 
nearly-elastic collisions, see Eq. (1), the dimensionless 
parameters e\ <~ e 2 ~ 1 — a 2 <C 1. The rescaled energy 
equation is 

dT 1 

— + (v • V)T = -(7 - 1)TV ■ v + -V ■ (VTVT) 



- 2pT 



3/2 . 



e 2 ( 7 -i) Vr 

2p 



S 2 . (A2) 



In our simulation we assumed a 2D granular gas of hard 
disks, so 7 = 2. 

The transformation we used is inspired by Haff 's law. 
We introduce a transformed time, temperature and ve- 
locity: 



r = ln(l + t), 
f = {l + t) 2 T, 
v= (l + t)v, 



(A3) 



whereas the (rescaled) density p remains the same. The 
transformed equations become 



dp 

8t 



+ V • (pv) = 



(A4) 
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— + (v-v)f = -(7-i)fv-v + -v-(Vr ; vf) 

OT p 

- 2pf^ + 2f+ e2il "^S. (A6) 

As one can see, two new source terms appear: the e\pv 
term in the transformed momentum equation, and the 2T 
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term in the transformed energy equation. Neither of the 
source terms involves any derivatives, and so they could 
be easily implemented in our numerical codes. Impor- 
tantly, a time-dependent cooling state with a stationary 
density in the original formulation corresponds to a true 
steady state, for all variables, in the transformed equa- 
tions. 

Two different hydrodynamic codes, Vulcan/ID (V1D) 
and Vulcan/2D (V2D), were used for a numerical solu- 
tion of Eqs. (A4)-(A6). The codes involve ALE (Arbi- 
trary Lagrange-Euler) schemes on staggered grids and 
have both explicit and implicit solvers for hydrodynam- 
ics and heat conduction [51]. V1D mainly works in the 
Lagrangian mode and therefore can efficiently handle ID 
density singularities. V2D is naturally limited to Eule- 
rian simulations. It can employ, however, a strongly non- 
uniform polar mesh. This feature of the code was very 
important when following the developing singularity in a 
circular container, see section IV B. The two codes were 
previously tested against other codes and used for a va- 
riety of astrophysical problems, as well as for granular 
hydrodynamic simulations [25-30]. 



of the gas varies as I free ~ 1/p ~ {t* — t) 2 . Therefore, 



lf r 



(t. - if 12 



0. 



(B2) 



and the applicability of hydrodynamics improves near the 
singularity. 

Finally, we will check whether the pressure profile re- 
mains locally uniform near the singularity. This amounts 
to checking whether the correction p^\r,t), which we 
neglected in the perturbative expansion (58), is indeed 
much smaller than the term p(°\t). To evaluate we 
return to Eq. (3) and plug in it the (supposedly leading- 
order) quantities found in section VI B: 



dr 



dv dv 



dt 



dr 




(B3) 



Let us estimate the terms on the right hand side. The 
viscous term behaves, in the vicinity of the vertex of the 
container r — 0, as 



Appendix B. Verifying assumptions a posteriori 

As we have seen in subsection VIC, heat diffusion is 
irrelevant for thermal collapse in 2D or 3D. Let us check 
other assumptions that we made in order to derive the 
scaling relations for thermal collapse. 

When using ideal equations, we neglected the viscous 
heating term in the energy equation (5) or (6). The ratio 
of this term to the energy loss term can be estimated as 



^0(7 



2Ap 1 / 2 p 3 / 2 



z/oA / p_ p_ 
2 V Po Ps 
1 - a 2 < 1 . 



(Bl) 



Therefore, for nearly elastic collisions (that we assume 
throughout this work) , the viscous heating remains small 
and, in the leading order, can be neglected. 

The validity of hydrodynamics itself at the density 
blowup demands that the Knudsen number remains 
small, see Eq. (1). As we have seen in section VI B, 
the smallest hydrodynamic scale varies with time, in 2D, 
as I ~ (t* — i) 3 / 2 . On the other hand, the mean free path 



^0 v .{^m/'pt)^ v j r ^Tp v 1 

VQ^Ps n 2,Ps „ Ps 

~ ■*T7~ (1 - a) 7 < 7- 



(B4) 



As d r p^ <~ p^ /£, the viscous term does not violate the 
strong inequality p^ -C p(°h 

Finally, let us estimate the inertial terms: 



f dv dv\ pv 2 p s (t* — t 



2/3-4 



• (B5) 



Therefore, the contribution of the inertial terms to p^ 
is such that 



P {1) ,, 2^ (U-l\ r 



(B6) 



One can see that, in the limit of nearly elastic collisions, 
1 — a 2 <C 1, the ratio p^/pW is small at early times. 
However, for d = 2 or 3 it grows without limit as t — > i». 
This suggests that the uniform-pressure approximation 
breaks down sufficiently close to the time of singularity. 
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